/*
 * BeastGenerator.java
 *
 * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut
 *
 * This file is part of BEAST.
 * See the NOTICE file distributed with this work for additional
 * information regarding copyright ownership and licensing.
 *
 * BEAST is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as
 * published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 *  BEAST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with BEAST; if not, write to the
 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
 * Boston, MA  02110-1301  USA
 */

package dr.app.beauti;

import dr.evolution.alignment.SitePatterns;
import dr.evolution.datatype.Nucleotides;
import dr.evolution.tree.NodeRef;
import dr.evolution.tree.Tree;
import dr.evolution.util.Taxon;
import dr.evolution.util.TaxonList;
import dr.evolution.util.Units;
import dr.evomodel.branchratemodel.DiscretizedBranchRates;
import dr.evomodel.branchratemodel.StrictClockBranchRates;
import dr.evomodel.coalescent.*;
import dr.evomodel.operators.ExchangeOperator;
import dr.evomodel.operators.SubtreeSlideOperator;
import dr.evomodel.operators.WilsonBalding;
import dr.evomodel.sitemodel.GammaSiteModel;
import dr.evomodel.speciation.BirthDeathModel;
import dr.evomodel.speciation.SpeciationLikelihood;
import dr.evomodel.speciation.YuleModel;
import dr.evomodel.substmodel.EmpiricalAminoAcidModel;
import dr.evomodel.substmodel.FrequencyModel;
import dr.evomodel.tree.RateCovarianceStatistic;
import dr.evomodel.tree.RateStatistic;
import dr.evomodel.tree.TreeLogger;
import dr.evomodel.tree.TreeModel;
import dr.evomodel.treelikelihood.TreeLikelihood;
import dr.evoxml.*;
import dr.inference.distribution.DistributionLikelihood;
import dr.inference.distribution.ExponentialMarkovModel;
import dr.inference.distribution.LogNormalDistributionModel;
import dr.inference.loggers.Columns;
import dr.inference.loggers.MCLogger;
import dr.inference.model.*;
import dr.inference.operators.*;
import dr.util.Attribute;

import java.io.Writer;
import java.util.ArrayList;
import java.util.HashSet;
import java.util.Iterator;
import java.util.Set;

/**
 * This class holds all the data for the current BEAUti Document
 *
 * @author			Andrew Rambaut
 * @author			Alexei Drummond
 * @version			$Id: BeastGenerator.java,v 1.4 2006/09/05 13:29:34 rambaut Exp $
 */
public class BeastGenerator extends BeautiOptions {

    public BeastGenerator() {
        super();
    }

    /**
     * Checks various options to check they are valid. Throws IllegalArgumentExceptions with
     * descriptions of the problems.
     * @throws IllegalArgumentException if there is a problem with the current settings
     */
    public void checkOptions() throws IllegalArgumentException {
        Set ids = new HashSet();

        ids.add("taxa");
        ids.add("alignment");

        if (taxonList != null) {
            for (int i = 0; i < taxonList.getTaxonCount(); i++) {
                Taxon taxon = (Taxon)taxonList.getTaxon(i);
                if (ids.contains(taxon.getId())) {
                    throw new IllegalArgumentException("A taxon has the same id," + taxon.getId() +
                            "\nas another element (taxon, sequence, taxon set etc.):\nAll ids should be unique.");
                }
                ids.add(taxon.getId());
            }
        }

        for (int i = 0; i < taxonSets.size(); i++) {
            TaxonList taxa = (TaxonList)taxonSets.get(i);
            if (taxa.getTaxonCount() < 2) {
                throw new IllegalArgumentException("Taxon set, " + taxa.getId() + ", should contain\n"+
                        "at least two taxa.");
            }
            if (ids.contains(taxa.getId())) {
                throw new IllegalArgumentException("A taxon sets has the same id," + taxa.getId() +
                        "\nas another element (taxon, sequence, taxon set etc.):\nAll ids should be unique.");
            }
            ids.add(taxa.getId());
        }

        getPartionCount(codonHeteroPattern);
    }

    /**
     * Generate a beast xml file from these beast options
     * @param w the writer
     */
    public void generateXML(Writer w) {

        XMLWriter writer = new XMLWriter(w);

        writer.writeText("<?xml version=\"1.0\" standalone=\"yes\"?>");
        writer.writeComment("Generated by BEAUTi v1.4");
        writer.writeComment("      by Alexei J. Drummond and Andrew Rambaut");
        writer.writeComment("      Department of Computer Science, University of Auckland and");
        writer.writeComment("      Institute of Evolutionary Biology, University of Edinburgh");
        writer.writeComment("      http://evolve.zoo.ox.ac.uk/beast/");
        writer.writeOpenTag("beast");
        writer.writeText("");
        writeTaxa(writer);

        if (taxonSets != null && taxonSets.size() > 0) {
            writeTaxonSets(writer);
        }

        if (alignment != null) {
            writeAlignment(writer);
            writePatternLists(writer);
        }

        writer.writeText("");
        writeNodeHeightPriorModel(writer);

        writer.writeText("");
        writeStartingTree(writer);
        writer.writeText("");
        writeTreeModel(writer);
        writer.writeText("");
        writeNodeHeightPrior(writer);
        if (nodeHeightPrior == LOGISTIC) {
            writer.writeText("");
            writeBooleanLikelihood(writer);
        } else if (nodeHeightPrior == SKYLINE) {
            writer.writeText("");
            writeExponentialMarkovLikelihood(writer);
        }

        writer.writeText("");
        writeBranchRatesModel(writer);

        if (alignment != null) {
            writer.writeText("");
            writeSubstitutionModel(writer);
            writer.writeText("");
            writeSiteModel(writer);
            writer.writeText("");
            writeTreeLikelihood(writer);
        }

        writer.writeText("");

        if (taxonSets != null && taxonSets.size() > 0) {
            writeTMRCAStatistics(writer);
        }

        ArrayList operators = selectOperators();
        writeOperatorSchedule(operators, writer);
        writer.writeText("");
        writeMCMC(writer);
        writer.writeText("");
        writeTimerReport(writer);
        writer.writeText("");
        if (performTraceAnalysis) {
            writeTraceAnalysis(writer);
        }

        writer.writeCloseTag("beast");
        writer.flush();
    }

    /**
     * Generate a taxa block from these beast options
     * @param writer the writer
     */
    public void writeTaxa(XMLWriter writer) {

        writer.writeComment("The list of taxa analyse (can also include dates/ages).");
        writer.writeComment("ntax=" + taxonList.getTaxonCount());
        writer.writeOpenTag("taxa", new Attribute[] { new Attribute.Default("id", "taxa") });

        boolean firstDate = true;
        for (int i = 0; i < taxonList.getTaxonCount(); i++) {
            Taxon taxon = taxonList.getTaxon(i);

            boolean hasDate = false;

            if (maximumTipHeight > 0.0) {
                hasDate = TaxonList.Utils.hasAttribute(taxonList, i, dr.evolution.util.Date.DATE);
            }

            writer.writeTag("taxon", new Attribute[] { new Attribute.Default("id", taxon.getId()) }, !hasDate);

            if (hasDate) {
                dr.evolution.util.Date date = (dr.evolution.util.Date)taxon.getAttribute(dr.evolution.util.Date.DATE);

                if (firstDate) {
                    units = date.getUnits();
                    firstDate = false;
                } else {
                    if (units != date.getUnits()) {
                        System.err.println("Error: Units in dates do not match.");
                    }
                }

                Attribute[] attributes = new Attribute[] {
                        new Attribute.Default("value", date.getTimeValue()+""),
                        new Attribute.Default("direction", date.isBackwards() ? "backwards" : "forwards"),
                        new Attribute.Default("units", Units.UNIT_NAMES[units][0])
                        /*,
                                                  new Attribute.Default("origin", date.getOrigin()+"")*/
                };

                writer.writeTag(dr.evolution.util.Date.DATE, attributes, true);
                writer.writeCloseTag("taxon");
            }
        }

        writer.writeCloseTag("taxa");
    }

    /**
     * Generate additional taxon sets
     * @param writer the writer
     */
    public void writeTaxonSets(XMLWriter writer) {

        writer.writeText("");
        for (int i = 0; i < taxonSets.size(); i++) {
            TaxonList taxa = (TaxonList)taxonSets.get(i);
            writer.writeOpenTag(
                    "taxa",
                    new Attribute[] {
                            new Attribute.Default("id", taxa.getId())
                    }
            );

            for (int j = 0; j < taxa.getTaxonCount(); j++) {
                Taxon taxon = taxa.getTaxon(j);

                writer.writeTag("taxon", new Attribute[] { new Attribute.Default("idref", taxon.getId()) }, true);
            }
            writer.writeCloseTag("taxa");
        }
    }

    /**
     * Generate an alignment block from these beast options
     * @param writer the writer
     */
    public void writeAlignment(XMLWriter writer) {

        writer.writeText("");
        writer.writeComment("The sequence alignment (each sequence refers to a taxon above).");
        writer.writeComment("ntax=" + alignment.getTaxonCount() + " nchar=" + alignment.getSiteCount());
        writer.writeOpenTag(
                "alignment",
                new Attribute[] {
                        new Attribute.Default("id", "alignment"),
                        new Attribute.Default("dataType", alignment.getDataType().getDescription())
                }
        );

        for (int i = 0; i < alignment.getTaxonCount(); i++) {
            Taxon taxon = alignment.getTaxon(i);

            writer.writeOpenTag("sequence");
            writer.writeTag("taxon", new Attribute[] { new Attribute.Default("idref", taxon.getId()) }, true);
            writer.writeText(alignment.getAlignedSequenceString(i));
            writer.writeCloseTag("sequence");
        }
        writer.writeCloseTag("alignment");
    }

    /**
     * Write a demographic model
     * @param writer the writer
     */
    public void writeNodeHeightPriorModel(XMLWriter writer) {

        String initialPopSize = null;

        if (nodeHeightPrior == CONSTANT) {

            writer.writeComment("A prior assumption that the population size has remained constant");
            writer.writeComment("throughout the time spanned by the genealogy.");
            writer.writeOpenTag(
                    ConstantPopulationModel.CONSTANT_POPULATION_MODEL,
                    new Attribute[] {
                            new Attribute.Default("id", "constant"),
                            new Attribute.Default("units", Units.Utils.getDefaultUnitName(units))
                    }
            );

            writer.writeOpenTag(ConstantPopulationModel.POPULATION_SIZE);
            writeParameter("constant.popSize", writer);
            writer.writeCloseTag(ConstantPopulationModel.POPULATION_SIZE);
            writer.writeCloseTag(ConstantPopulationModel.CONSTANT_POPULATION_MODEL);

        } else if (nodeHeightPrior == EXPONENTIAL) {
            // generate an exponential prior tree

            writer.writeComment("A prior assumption that the population size has grown exponentially");
            writer.writeComment("throughout the time spanned by the genealogy.");
            writer.writeOpenTag(
                    ExponentialGrowthModel.EXPONENTIAL_GROWTH_MODEL,
                    new Attribute[] {
                            new Attribute.Default("id", "exponential"),
                            new Attribute.Default("units", Units.Utils.getDefaultUnitName(units))
                    }
            );

            // write pop size socket
            writer.writeOpenTag(ExponentialGrowthModel.POPULATION_SIZE);
            writeParameter("exponential.popSize", writer);
            writer.writeCloseTag(ExponentialGrowthModel.POPULATION_SIZE);

            if (parameterization == GROWTH_RATE) {
                // write growth rate socket
                writer.writeOpenTag(ExponentialGrowthModel.GROWTH_RATE);
                writeParameter("exponential.growthRate", writer);
                writer.writeCloseTag(ExponentialGrowthModel.GROWTH_RATE);
            } else {
                // write doubling time socket
                writer.writeOpenTag(ExponentialGrowthModel.DOUBLING_TIME);
                writeParameter("exponential.doublingTime", writer);
                writer.writeCloseTag(ExponentialGrowthModel.DOUBLING_TIME);
            }

            writer.writeCloseTag(ExponentialGrowthModel.EXPONENTIAL_GROWTH_MODEL);
        } else if (nodeHeightPrior == LOGISTIC) {
            // generate an exponential prior tree

            writer.writeComment("A prior assumption that the population size has grown logistically");
            writer.writeComment("throughout the time spanned by the genealogy.");
            writer.writeOpenTag(
                    LogisticGrowthModel.LOGISTIC_GROWTH_MODEL,
                    new Attribute[] {
                            new Attribute.Default("id", "logistic"),
                            new Attribute.Default("units", Units.Utils.getDefaultUnitName(units))
                    }
            );

            // write pop size socket
            writer.writeOpenTag(LogisticGrowthModel.POPULATION_SIZE);
            writeParameter("logistic.popSize", writer);
            writer.writeCloseTag(LogisticGrowthModel.POPULATION_SIZE);

            if (parameterization == GROWTH_RATE) {
                // write growth rate socket
                writer.writeOpenTag(LogisticGrowthModel.GROWTH_RATE);
                writeParameter("logistic.growthRate", writer);
                writer.writeCloseTag(LogisticGrowthModel.GROWTH_RATE);
            } else {
                // write doubling time socket
                writer.writeOpenTag(LogisticGrowthModel.DOUBLING_TIME);
                writeParameter("logistic.doublingTime", writer);
                writer.writeCloseTag(LogisticGrowthModel.DOUBLING_TIME);
            }

            // write logistic t50 socket
            writer.writeOpenTag(LogisticGrowthModel.TIME_50);
            writeParameter("logistic.t50", writer);
            writer.writeCloseTag(LogisticGrowthModel.TIME_50);

            writer.writeCloseTag(LogisticGrowthModel.LOGISTIC_GROWTH_MODEL);

            initialPopSize = "logistic.popSize";

        } else if (nodeHeightPrior == EXPANSION) {
            // generate an exponential prior tree

            writer.writeComment("A prior assumption that the population size has grown exponentially");
            writer.writeComment("from some ancestral population size in the past.");
            writer.writeOpenTag(
                    ExpansionModel.EXPANSION_MODEL,
                    new Attribute[] {
                            new Attribute.Default("id", "expansion"),
                            new Attribute.Default("units", Units.Utils.getDefaultUnitName(units))
                    }
            );

            // write pop size socket
            writer.writeOpenTag(ExpansionModel.POPULATION_SIZE);
            writeParameter("expansion.popSize", writer);
            writer.writeCloseTag(ExpansionModel.POPULATION_SIZE);

            if (parameterization == GROWTH_RATE) {
                // write growth rate socket
                writer.writeOpenTag(ExpansionModel.GROWTH_RATE);
                writeParameter("expansion.growthRate", writer);
                writer.writeCloseTag(ExpansionModel.GROWTH_RATE);
            } else {
                // write doubling time socket
                writer.writeOpenTag(ExpansionModel.DOUBLING_TIME);
                writeParameter("expansion.doublingTime", writer);
                writer.writeCloseTag(ExpansionModel.DOUBLING_TIME);
            }

            // write ancestral proportion socket
            writer.writeOpenTag(ExpansionModel.ANCESTRAL_POPULATION_PROPORTION);
            writeParameter("expansion.ancestralProportion", writer);
            writer.writeCloseTag(ExpansionModel.ANCESTRAL_POPULATION_PROPORTION);

            writer.writeCloseTag(ExpansionModel.EXPANSION_MODEL);

            initialPopSize = "expansion.popSize";

        } else if (nodeHeightPrior == YULE) {
            writer.writeText("");
            writer.writeComment("A prior on the distribution node heights defined given");
            writer.writeComment("a Yule speciation process (a pure birth process).");
            writer.writeOpenTag(
                    YuleModel.YULE_MODEL,
                    new Attribute[] {
                            new Attribute.Default("id", "yule"),
                            new Attribute.Default("units", Units.Utils.getDefaultUnitName(units))
                    }
            );

            writer.writeOpenTag(YuleModel.BIRTH_RATE);
            writeParameter("yule.birthRate", writer);
            writer.writeCloseTag(YuleModel.BIRTH_RATE);
            writer.writeCloseTag(YuleModel.YULE_MODEL);
        } else if (nodeHeightPrior == BIRTH_DEATH) {
            writer.writeText("");
            writer.writeComment("A prior on the distribution node heights defined given");
            writer.writeComment("a Birth-Death speciation process (Yang & Rannala, 1997).");
            writer.writeOpenTag(
                    BirthDeathModel.BIRTH_DEATH_MODEL,
                    new Attribute[] {
                            new Attribute.Default("id", "birthDeath"),
                            new Attribute.Default("units", Units.Utils.getDefaultUnitName(units))
                    }
            );

            writer.writeOpenTag(BirthDeathModel.BIRTH_RATE);
            writeParameter("birthDeath.birthRate", writer);
            writer.writeCloseTag(BirthDeathModel.BIRTH_RATE);
            writer.writeOpenTag(BirthDeathModel.DEATH_RATE);
            writeParameter("birthDeath.deathRate", writer);
            writer.writeCloseTag(BirthDeathModel.DEATH_RATE);
            writer.writeOpenTag(BirthDeathModel.SAMPLING_PROPORTION);
            // We are not sampling this parameter so use a fixed value
            writeParameter("birthDeath.samplingProportion", 1, birthDeathSamplingProportion, Double.NaN, Double.NaN, writer);
            writer.writeCloseTag(BirthDeathModel.SAMPLING_PROPORTION);
            writer.writeCloseTag(BirthDeathModel.BIRTH_DEATH_MODEL);
        }

        if (nodeHeightPrior != CONSTANT && nodeHeightPrior != EXPONENTIAL) {
            // If the node height prior is not one of these two then we need to simulate a
            // random starting tree under a constant size coalescent.

            writer.writeComment("This is a simple constant population size coalescent model");
            writer.writeComment("that is used to generate an initial tree for the chain.");
            writer.writeOpenTag(
                    ConstantPopulationModel.CONSTANT_POPULATION_MODEL,
                    new Attribute[] {
                            new Attribute.Default("id", "initialDemo"),
                            new Attribute.Default("units", Units.Utils.getDefaultUnitName(units))
                    }
            );

            writer.writeOpenTag(ConstantPopulationModel.POPULATION_SIZE);
            if (initialPopSize != null) {
                writer.writeTag(ParameterParser.PARAMETER,
                        new Attribute[] {
                                new Attribute.Default("idref", initialPopSize),
                        }, true);
            } else {
                writeParameter("initialDemo.popSize", 1, 100.0, Double.NaN, Double.NaN, writer);
            }
            writer.writeCloseTag(ConstantPopulationModel.POPULATION_SIZE);
            writer.writeCloseTag(ConstantPopulationModel.CONSTANT_POPULATION_MODEL);
        }

    }

    /**
     * Writes the pattern lists
     * @param writer the writer
     */
    public void writePatternLists(XMLWriter writer) {

        partitionCount = getPartionCount(codonHeteroPattern);

        writer.writeText("");
        if (alignment.getDataType() == Nucleotides.INSTANCE && codonHeteroPattern != null && partitionCount > 1) {

            if (codonHeteroPattern.equals("112")) {
                writer.writeComment("The unique patterns for codon positions 1 & 2");
                writer.writeOpenTag(MergePatternsParser.MERGE_PATTERNS,
                        new Attribute[] {
                                new Attribute.Default("id", "patterns1+2"),
                        }
                );
                writePatternList(1, 3, writer);
                writePatternList(2, 3, writer);
                writer.writeCloseTag(MergePatternsParser.MERGE_PATTERNS);

                writePatternList(3, 3, writer);

            } else {
                // pattern is 123
                // write pattern lists for all three codon positions
                for (int i = 1; i <= 3; i++) {
                    writePatternList(i, 3, writer);
                }

            }
        } else {
            partitionCount = 1;
            writePatternList(-1, 0, writer);
        }
    }

    private int getPartionCount(String codonPattern) {

        if (codonPattern == null || codonPattern.equals("111")) {
            return 1;
        }
        if (codonPattern.equals("123")) {
            return 3;
        }
        if (codonPattern.equals("112")) {
            return 2;
        }
        throw new IllegalArgumentException("codonPattern must be one of '111', '112' or '123'");
    }

    /**
     * Write a single pattern list
     * @param writer the writer
     * @param from from site
     * @param every skip every
     */
    private void writePatternList(int from, int every, XMLWriter writer) {

        String id = "patterns";
        if (from < 1) {
            writer.writeComment("The unique patterns for all positions");
            from = 1;
        } else {
            writer.writeComment("The unique patterns for codon position "+from);
            id += Integer.toString(from);
        }

        SitePatterns patterns = new SitePatterns(alignment, from-1, 0, every);
        writer.writeComment("npatterns=" + patterns.getPatternCount());
        if (every != 0) {
            writer.writeOpenTag(SitePatternsParser.PATTERNS,
                    new Attribute[] {
                            new Attribute.Default("id", id),
                            new Attribute.Default("from", ""+from),
                            new Attribute.Default("every", ""+every)
                    }
            );
        } else {
            writer.writeOpenTag(SitePatternsParser.PATTERNS,
                    new Attribute[] {
                            new Attribute.Default("id", id),
                            new Attribute.Default("from", ""+from)
                    }
            );
        }

        writer.writeTag("alignment", new Attribute.Default("idref", "alignment"), true);
        writer.writeCloseTag(SitePatternsParser.PATTERNS);
    }

    /**
     * Write tree model XML block.
     * @param writer the writer
     */
    private void writeTreeModel(XMLWriter writer) {

        writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("id", "treeModel"), false);

        if (userTree) {
            writer.writeTag("tree", new Attribute.Default("idref", "startingTree"), true);
        } else {
            writer.writeTag(CoalescentSimulator.COALESCENT_TREE, new Attribute.Default("idref", "startingTree"), true);
        }

        writer.writeOpenTag(TreeModel.ROOT_HEIGHT);
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("id", "treeModel.rootHeight"), true);
        writer.writeCloseTag(TreeModel.ROOT_HEIGHT);


        writer.writeOpenTag(TreeModel.NODE_HEIGHTS, new Attribute.Default(TreeModel.INTERNAL_NODES, "true"));
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("id", "treeModel.internalNodeHeights"), true);
        writer.writeCloseTag(TreeModel.NODE_HEIGHTS);

        writer.writeOpenTag(TreeModel.NODE_HEIGHTS,
                new Attribute[] {
                        new Attribute.Default(TreeModel.INTERNAL_NODES, "true"),
                        new Attribute.Default(TreeModel.ROOT_NODE, "true")
                });
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("id", "treeModel.allInternalNodeHeights"), true);
        writer.writeCloseTag(TreeModel.NODE_HEIGHTS);

        writer.writeCloseTag(TreeModel.TREE_MODEL);
    }

    /**
     * Writes the substitution model to XML.
     * @param writer the writer
     */
    public void writeSubstitutionModel(XMLWriter writer) {

        if (alignment.getDataType() == Nucleotides.INSTANCE) {

            // Jukes-Cantor model
            if (nucSubstitutionModel == JC) {
                writer.writeComment("The JC substitution model (Jukes & Cantor, 1969)");
                writer.writeOpenTag(
                        dr.evomodel.substmodel.HKY.HKY_MODEL,
                        new Attribute[] { new Attribute.Default("id", "jc") }
                );
                writer.writeOpenTag(dr.evomodel.substmodel.HKY.FREQUENCIES);
                writer.writeOpenTag(
                        FrequencyModel.FREQUENCY_MODEL,
                        new Attribute[] {
                                new Attribute.Default("dataType", alignment.getDataType().getDescription())
                        }
                );
                writer.writeOpenTag(FrequencyModel.FREQUENCIES);
                writer.writeTag(
                        ParameterParser.PARAMETER,
                        new Attribute[] {
                                new Attribute.Default("id", "jc.frequencies"),
                                new Attribute.Default("value", "0.25 0.25 0.25 0.25")
                        },
                        true
                );
                writer.writeCloseTag(FrequencyModel.FREQUENCIES);

                writer.writeCloseTag(FrequencyModel.FREQUENCY_MODEL);
                writer.writeCloseTag(dr.evomodel.substmodel.HKY.FREQUENCIES);

                writer.writeOpenTag(dr.evomodel.substmodel.HKY.KAPPA);
                writeParameter("jc.kappa", 1, 1.0, Double.NaN, Double.NaN, writer);
                writer.writeCloseTag(dr.evomodel.substmodel.HKY.KAPPA);
                writer.writeCloseTag(dr.evomodel.substmodel.HKY.HKY_MODEL);

            } else {
                // Hasegawa Kishino and Yano 85 model
                if (nucSubstitutionModel == HKY) {
                    if (unlinkedSubstitutionModel) {
                        for (int i = 1; i <= partitionCount; i++) {
                            writeHKYModel(i, writer);
                        }
                    } else {
                        writeHKYModel(-1, writer);
                    }
                } else {
                    // General time reversible model
                    if (nucSubstitutionModel == GTR) {
                        if (unlinkedSubstitutionModel) {
                            for (int i = 1; i <= partitionCount; i++) {
                                writeGTRModel(i, writer);
                            }
                        } else {
                            writeGTRModel(-1, writer);
                        }
                    }
                }
            }
        } else {
            // Amino Acid model
            String aaModel = "";

            switch (aaSubstitutionModel) {
                case 0: aaModel = EmpiricalAminoAcidModel.BLOSUM_62; break;
                case 1: aaModel = EmpiricalAminoAcidModel.DAYHOFF; break;
                case 2: aaModel = EmpiricalAminoAcidModel.JTT; break;
                case 3: aaModel = EmpiricalAminoAcidModel.MT_REV_24; break;
                case 4: aaModel = EmpiricalAminoAcidModel.CP_REV_45; break;
                case 5: aaModel = EmpiricalAminoAcidModel.WAG; break;
            }

            writer.writeComment("The " + aaModel + " substitution model");
            writer.writeTag(
                    EmpiricalAminoAcidModel.EMPIRICAL_AMINO_ACID_MODEL,
                    new Attribute[] { new Attribute.Default("id", "aa"),
                            new Attribute.Default("type", aaModel) }, true
            );

        }
    }

    /**
     * Write the HKY model XML block.
     * @param num the model number
     * @param writer the writer
     */
    public void writeHKYModel(int num, XMLWriter writer) {
        String id = "hky";
        if (num > 0) {
            id += Integer.toString(num);
        }
        // Hasegawa Kishino and Yano 85 model
        writer.writeComment("The HKY substitution model (Hasegawa, Kishino & Yano, 1985)");
        writer.writeOpenTag(
                dr.evomodel.substmodel.HKY.HKY_MODEL,
                new Attribute[] { new Attribute.Default("id", id) }
        );
        writer.writeOpenTag(dr.evomodel.substmodel.HKY.FREQUENCIES);
        writer.writeOpenTag(
                FrequencyModel.FREQUENCY_MODEL,
                new Attribute[] {
                        new Attribute.Default("dataType", alignment.getDataType().getDescription())
                }
        );
        writer.writeTag("alignment", new Attribute[] { new Attribute.Default("idref", "alignment") }, true);
        writer.writeOpenTag(FrequencyModel.FREQUENCIES);
        writer.writeTag(
                ParameterParser.PARAMETER,
                new Attribute[] {
                        new Attribute.Default("id", id + ".frequencies"),
                        new Attribute.Default("dimension", "4")
                },
                true
        );
        writer.writeCloseTag(FrequencyModel.FREQUENCIES);
        writer.writeCloseTag(FrequencyModel.FREQUENCY_MODEL);
        writer.writeCloseTag(dr.evomodel.substmodel.HKY.FREQUENCIES);

        writer.writeOpenTag(dr.evomodel.substmodel.HKY.KAPPA);
        writeParameter(id + ".kappa", writer);
        writer.writeCloseTag(dr.evomodel.substmodel.HKY.KAPPA);
        writer.writeCloseTag(dr.evomodel.substmodel.HKY.HKY_MODEL);
    }

    /**
     * Write the GTR model XML block.
     * @param num the model number
     * @param writer the writer
     */
    public void writeGTRModel(int num, XMLWriter writer) {
        String id = "gtr";
        if (num > 0) {
            id += Integer.toString(num);
        }

        writer.writeComment("The general time reversible (GTR) substitution model");
        writer.writeOpenTag(
                dr.evomodel.substmodel.GTR.GTR_MODEL,
                new Attribute[] { new Attribute.Default("id", id) }
        );
        writer.writeOpenTag(dr.evomodel.substmodel.GTR.FREQUENCIES);
        writer.writeOpenTag(
                FrequencyModel.FREQUENCY_MODEL,
                new Attribute[] {
                        new Attribute.Default("dataType", alignment.getDataType().getDescription())
                }
        );
        writer.writeTag("alignment", new Attribute[] { new Attribute.Default("idref", "alignment") }, true);
        writer.writeOpenTag(FrequencyModel.FREQUENCIES);
        writer.writeTag(
                ParameterParser.PARAMETER,
                new Attribute[] {
                        new Attribute.Default("id", id + ".frequencies"),
                        new Attribute.Default("dimension", "4")
                },
                true
        );
        writer.writeCloseTag(FrequencyModel.FREQUENCIES);
        writer.writeCloseTag(FrequencyModel.FREQUENCY_MODEL);
        writer.writeCloseTag(dr.evomodel.substmodel.GTR.FREQUENCIES);

        writer.writeOpenTag(dr.evomodel.substmodel.GTR.A_TO_C);
        writeParameter(id + ".ac", writer);
        writer.writeCloseTag(dr.evomodel.substmodel.GTR.A_TO_C);

        writer.writeOpenTag(dr.evomodel.substmodel.GTR.A_TO_G);
        writeParameter(id + ".ag", writer);
        writer.writeCloseTag(dr.evomodel.substmodel.GTR.A_TO_G);

        writer.writeOpenTag(dr.evomodel.substmodel.GTR.A_TO_T);
        writeParameter(id + ".at", writer);
        writer.writeCloseTag(dr.evomodel.substmodel.GTR.A_TO_T);

        writer.writeOpenTag(dr.evomodel.substmodel.GTR.C_TO_G);
        writeParameter(id + ".cg", writer);
        writer.writeCloseTag(dr.evomodel.substmodel.GTR.C_TO_G);

        writer.writeOpenTag(dr.evomodel.substmodel.GTR.G_TO_T);
        writeParameter(id + ".gt", writer);
        writer.writeCloseTag(dr.evomodel.substmodel.GTR.G_TO_T);
        writer.writeCloseTag(dr.evomodel.substmodel.GTR.GTR_MODEL);
    }

    /**
     * Write the site model XML block.
     * @param writer the writer
     */
    public void writeSiteModel(XMLWriter writer) {
        if (alignment.getDataType() == Nucleotides.INSTANCE) {
            if (codonHeteroPattern != null) {
                for (int i = 1; i <= partitionCount; i++) {
                    writeNucSiteModel(i, writer);
                }
                writer.println();
                writer.writeOpenTag(CompoundParameter.COMPOUND_PARAMETER, new Attribute[] {new Attribute.Default("id", "allMus")});
                for (int i = 1; i <= partitionCount; i++) {
                    writer.writeTag(ParameterParser.PARAMETER,
                            new Attribute[] {new Attribute.Default("idref", "siteModel" + i + ".mu")}, true);
                }
                writer.writeCloseTag(CompoundParameter.COMPOUND_PARAMETER);
            } else {
                writeNucSiteModel(-1, writer);
            }
        } else {
            writeAASiteModel(writer);
        }
    }

    /**
     * Write the nucleotide site model XML block.
     * @param num the model number
     * @param writer the writer
     */
    public void writeNucSiteModel(int num, XMLWriter writer) {

        String id = "siteModel";
        if (num > 0) {
            id += Integer.toString(num);
        }

        writer.writeComment("site model");
        writer.writeOpenTag(GammaSiteModel.SITE_MODEL,new Attribute[] { new Attribute.Default("id", id) });


        writer.writeOpenTag(GammaSiteModel.SUBSTITUTION_MODEL);

        if (unlinkedSubstitutionModel) {
            switch (nucSubstitutionModel) {
                // JC cannot be unlinked because it has no parameters
                case JC: writer.writeTag(dr.evomodel.substmodel.HKY.HKY_MODEL, new Attribute.Default("idref", "jc"), true); break;
                case HKY: writer.writeTag(dr.evomodel.substmodel.HKY.HKY_MODEL, new Attribute.Default("idref", "hky" + num), true); break;
                case GTR: writer.writeTag(dr.evomodel.substmodel.GTR.GTR_MODEL, new Attribute.Default("idref", "gtr" + num), true); break;
                default: throw new IllegalArgumentException("Unknown substitution model.");
            }
        } else {
            switch (nucSubstitutionModel) {
                case JC: writer.writeTag(dr.evomodel.substmodel.HKY.HKY_MODEL, new Attribute.Default("idref", "jc"), true); break;
                case HKY: writer.writeTag(dr.evomodel.substmodel.HKY.HKY_MODEL, new Attribute.Default("idref", "hky"), true); break;
                case GTR: writer.writeTag(dr.evomodel.substmodel.GTR.GTR_MODEL, new Attribute.Default("idref", "gtr"), true); break;
                default: throw new IllegalArgumentException("Unknown substitution model.");
            }
        }
        writer.writeCloseTag(GammaSiteModel.SUBSTITUTION_MODEL);

        if (num != -1) {
            writer.writeOpenTag(GammaSiteModel.RELATIVE_RATE);
            writeParameter(id + ".mu", writer);
            writer.writeCloseTag(GammaSiteModel.RELATIVE_RATE);
        } else {
//            The actual mutation rate is now in the BranchRateModel so relativeRate can be missing
        }

        if (gammaHetero) {
            writer.writeOpenTag(GammaSiteModel.GAMMA_SHAPE, new Attribute.Default(GammaSiteModel.GAMMA_CATEGORIES, ""+gammaCategories));
            if (num == -1 || unlinkedHeterogeneityModel) {
                writeParameter(id + ".alpha", writer);
            }  else {
                // multiple partitions but linked heterogeneity
                if (num == 1) {
                    writeParameter("siteModel.alpha", writer);
                } else {
                    writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", "siteModel.alpha"), true);
                }
            }
            writer.writeCloseTag(GammaSiteModel.GAMMA_SHAPE);
        }

        if (invarHetero) {
            writer.writeOpenTag(GammaSiteModel.PROPORTION_INVARIANT);
            if (num == -1 || unlinkedHeterogeneityModel) {
                writeParameter(id + ".pInv", writer);
            }  else {
                // multiple partitions but linked heterogeneity
                if (num == 1) {
                    writeParameter("siteModel.pInv", writer);
                } else {
                    writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", "siteModel.pInv"), true);
                }
            }
            writer.writeCloseTag(GammaSiteModel.PROPORTION_INVARIANT);
        }

        writer.writeCloseTag(GammaSiteModel.SITE_MODEL);
    }

    /**
     * Write the AA site model XML block.
     * @param writer the writer
     */
    public void writeAASiteModel(XMLWriter writer) {

        writer.writeComment("site model");
        writer.writeOpenTag(GammaSiteModel.SITE_MODEL,new Attribute[] {
                new Attribute.Default("id", "siteModel") });


        writer.writeOpenTag(GammaSiteModel.SUBSTITUTION_MODEL);
        writer.writeTag(EmpiricalAminoAcidModel.EMPIRICAL_AMINO_ACID_MODEL,
                new Attribute.Default("idref", "aa"), true);
        writer.writeCloseTag(GammaSiteModel.SUBSTITUTION_MODEL);

//            The actual mutation rate is now in the BranchRateModel so relativeRate can be missing

        if (gammaHetero) {
            writer.writeOpenTag(GammaSiteModel.GAMMA_SHAPE, new Attribute.Default(GammaSiteModel.GAMMA_CATEGORIES, ""+gammaCategories));
            writeParameter("siteModel.alpha", writer);
            writer.writeCloseTag(GammaSiteModel.GAMMA_SHAPE);
        }

        if (invarHetero) {
            writer.writeOpenTag(GammaSiteModel.PROPORTION_INVARIANT);
            writeParameter("siteModel.pInv", writer);
            writer.writeCloseTag(GammaSiteModel.PROPORTION_INVARIANT);
        }

        writer.writeCloseTag(GammaSiteModel.SITE_MODEL);
    }

    /**
     * Write the relaxed clock branch rates block.
     * @param writer the writer
     */
    public void writeBranchRatesModel(XMLWriter writer) {
        if (clockModel == STRICT_CLOCK) {
            if (fixedSubstitutionRate) {
                fixParameter("clock.rate", meanSubstitutionRate);
            }

            writer.writeComment("The strict clock (Uniform rates across branches)");
            writer.writeOpenTag(
                    StrictClockBranchRates.STRICT_CLOCK_BRANCH_RATES,
                    new Attribute[] { new Attribute.Default("id", "branchRates") }
            );
            writer.writeOpenTag("rate");

            writeParameter("clock.rate", writer);
            writer.writeCloseTag("rate");
            writer.writeCloseTag(StrictClockBranchRates.STRICT_CLOCK_BRANCH_RATES);
        } else {
            writer.writeComment("The uncorrelated relaxed clock (Drummond, Ho, Phillips & Rambaut, 2005)");
            writer.writeOpenTag(
                    DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES,
                    new Attribute[] { new Attribute.Default("id", "branchRates") }
            );
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            writer.writeOpenTag("distribution");
            if (clockModel == UNCORRELATED_EXPONENTIAL) {
                if (fixedSubstitutionRate) {
                    fixParameter("uced.mean", meanSubstitutionRate);
                }

                writer.writeOpenTag("exponentialDistributionModel");
                writer.writeOpenTag("mean");
                writeParameter("uced.mean", writer);
                writer.writeCloseTag("mean");
                writer.writeCloseTag("exponentialDistributionModel");
            } else if (clockModel == UNCORRELATED_LOGNORMAL) {
                if (fixedSubstitutionRate) {
                    fixParameter("ucld.mean", meanSubstitutionRate);
                }

                writer.writeOpenTag("logNormalDistributionModel", new Attribute.Default(LogNormalDistributionModel.MEAN_IN_REAL_SPACE, "true"));
                writer.writeOpenTag("mean");
                writeParameter("ucld.mean", writer);
                writer.writeCloseTag("mean");
                writer.writeOpenTag("stdev");
                writeParameter("ucld.stdev", writer);
                writer.writeCloseTag("stdev");
                writer.writeCloseTag("logNormalDistributionModel");
            } else {
                throw new RuntimeException("Unrecognised relaxed clock model");
            }
            writer.writeCloseTag("distribution");
            writer.writeOpenTag("rateCategories");
            int categoryCount = (alignment.getSequenceCount() - 1) * 2;
            writeParameter("branchRates.categories", categoryCount, writer);
            writer.writeCloseTag("rateCategories");
            writer.writeCloseTag(DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES);

            writer.writeText("");
            writer.writeOpenTag(
                    RateStatistic.RATE_STATISTIC,
                    new Attribute[] {
                            new Attribute.Default("id", "meanRate"),
                            new Attribute.Default("name", "meanRate"),
                            new Attribute.Default("mode", "mean"),
                            new Attribute.Default("internal", "true"),
                            new Attribute.Default("external", "true")
                    }
            );
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            writer.writeTag(DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES, new Attribute.Default("idref", "branchRates"), true);
            writer.writeCloseTag(RateStatistic.RATE_STATISTIC);

            writer.writeText("");
            writer.writeOpenTag(
                    RateStatistic.RATE_STATISTIC,
                    new Attribute[] {
                            new Attribute.Default("id", "coefficientOfVariation"),
                            new Attribute.Default("name", "coefficientOfVariation"),
                            new Attribute.Default("mode", "coefficientOfVariation"),
                            new Attribute.Default("internal", "true"),
                            new Attribute.Default("external", "true")
                    }
            );
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            writer.writeTag(DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES, new Attribute.Default("idref", "branchRates"), true);
            writer.writeCloseTag(RateStatistic.RATE_STATISTIC);

            writer.writeText("");
            writer.writeOpenTag(
                    RateCovarianceStatistic.RATE_COVARIANCE_STATISTIC,
                    new Attribute[] {
                            new Attribute.Default("id", "covariance"),
                            new Attribute.Default("name", "covariance")
                    }
            );
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            writer.writeTag(DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES, new Attribute.Default("idref", "branchRates"), true);
            writer.writeCloseTag(RateCovarianceStatistic.RATE_COVARIANCE_STATISTIC);
        }
    }

    /**
     * Write the prior on node heights (coalescent or speciational models)
     * @param writer the writer
     */
    public void writeNodeHeightPrior(XMLWriter writer) {
        if (nodeHeightPrior == YULE || nodeHeightPrior == BIRTH_DEATH) {
            // generate a speciational process

            writer.writeOpenTag(
                    SpeciationLikelihood.SPECIATION_LIKELIHOOD,
                    new Attribute[] {
                            new Attribute.Default("id", "speciation")
                    }
            );

            // write pop size socket
            writer.writeOpenTag(SpeciationLikelihood.MODEL);
            writeNodeHeightPriorModelRef(writer);
            writer.writeCloseTag(SpeciationLikelihood.MODEL);
            writer.writeOpenTag(SpeciationLikelihood.TREE);
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            writer.writeCloseTag(SpeciationLikelihood.TREE);

            writer.writeCloseTag(SpeciationLikelihood.SPECIATION_LIKELIHOOD);

        } else if (nodeHeightPrior == SKYLINE) {
            // generate a Bayesian skyline plot

            writer.writeOpenTag(
                    BayesianSkylineLikelihood.SKYLINE_LIKELIHOOD,
                    new Attribute[] {
                            new Attribute.Default("id", "skyline"),
                            new Attribute.Default("linear", skylineModel == LINEAR_SKYLINE ? "true" : "false")
                    }
            );

            // write pop size socket
            writer.writeOpenTag(BayesianSkylineLikelihood.POPULATION_SIZES);
            if (skylineModel == LINEAR_SKYLINE) {
                writeParameter("skyline.popSize", skylineGroupCount + 1, writer);
            } else {
                writeParameter("skyline.popSize", skylineGroupCount, writer);
            }
            writer.writeCloseTag(BayesianSkylineLikelihood.POPULATION_SIZES);

            // write group size socket
            writer.writeOpenTag(BayesianSkylineLikelihood.GROUP_SIZES);
            writeParameter("skyline.groupSize", skylineGroupCount, writer);
            writer.writeCloseTag(BayesianSkylineLikelihood.GROUP_SIZES);

            writer.writeOpenTag(CoalescentLikelihood.POPULATION_TREE);
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            writer.writeCloseTag(CoalescentLikelihood.POPULATION_TREE);

            writer.writeCloseTag(BayesianSkylineLikelihood.SKYLINE_LIKELIHOOD);
        } else {
            // generate a coalescent process

            writer.writeOpenTag(
                    CoalescentLikelihood.COALESCENT_LIKELIHOOD,
                    new Attribute[] { new Attribute.Default("id", "coalescent") }
            );
            writer.writeOpenTag(CoalescentLikelihood.MODEL);
            writeNodeHeightPriorModelRef(writer);
            writer.writeCloseTag(CoalescentLikelihood.MODEL);
            writer.writeOpenTag(CoalescentLikelihood.POPULATION_TREE);
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            writer.writeCloseTag(CoalescentLikelihood.POPULATION_TREE);
            writer.writeCloseTag(CoalescentLikelihood.COALESCENT_LIKELIHOOD);
        }
    }

    /**
     * Write the boolean likelihood
     * @param writer the writer
     */
    public void writeBooleanLikelihood(XMLWriter writer) {
        writer.writeOpenTag(
                BooleanLikelihood.BOOLEAN_LIKELIHOOD,
                new Attribute[] { new Attribute.Default("id", "booleanLikelihood1") }
        );
        writer.writeOpenTag(
                TestStatistic.TEST_STATISTIC,
                new Attribute[] {
                        new Attribute.Default("id", "test1"),
                        new Attribute.Default("name", "test1")
                }
        );
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", "logistic.t50"), true);
        writer.writeOpenTag("lessThan");
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", "treeModel.rootHeight"), true);
        writer.writeCloseTag("lessThan");
        writer.writeCloseTag(TestStatistic.TEST_STATISTIC);
        writer.writeCloseTag(BooleanLikelihood.BOOLEAN_LIKELIHOOD);
    }

    public void writeExponentialMarkovLikelihood(XMLWriter writer) {
        writer.writeOpenTag(
                ExponentialMarkovModel.EXPONENTIAL_MARKOV_MODEL,
                new Attribute[] { new Attribute.Default("id", "eml1"),
                        new Attribute.Default("jeffreys", "true") }
        );
        writer.writeOpenTag( ExponentialMarkovModel.CHAIN_PARAMETER );
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", "skyline.popSize"), true);
        writer.writeCloseTag(ExponentialMarkovModel.CHAIN_PARAMETER);
        writer.writeCloseTag(ExponentialMarkovModel.EXPONENTIAL_MARKOV_MODEL);
    }

    /**
     * Write the tree likelihood XML block.
     * @param writer the writer
     */
    public void writeTreeLikelihood(XMLWriter writer) {
        boolean nucs = alignment.getDataType() == Nucleotides.INSTANCE;
        if (nucs && codonHeteroPattern != null) {
            for (int i = 1; i <= partitionCount; i++) {
                writeTreeLikelihood(i, writer);
            }
        } else {
            writeTreeLikelihood(-1, writer);
        }
    }

    /**
     * Write the tree likelihood XML block.
     * @param num the likelihood number
     * @param writer the writer
     */
    public void writeTreeLikelihood(int num, XMLWriter writer) {

        if (num > 0) {
            writer.writeOpenTag(
                    TreeLikelihood.TREE_LIKELIHOOD,
                    new Attribute[] { new Attribute.Default("id", "treeLikelihood" + num) }
            );
            if (codonHeteroPattern.equals("112") ) {
                if (num == 1) {
                    writer.writeTag(SitePatternsParser.PATTERNS, new Attribute[] { new Attribute.Default("idref", "patterns1+2") }, true);
                } else {
                    writer.writeTag(SitePatternsParser.PATTERNS, new Attribute[] { new Attribute.Default("idref", "patterns3") }, true);
                }
            } else {
                writer.writeTag(SitePatternsParser.PATTERNS, new Attribute[] { new Attribute.Default("idref", "patterns" + num) }, true);
            }
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute[] { new Attribute.Default("idref", "treeModel") }, true);
            writer.writeTag(GammaSiteModel.SITE_MODEL, new Attribute[] { new Attribute.Default("idref", "siteModel" + num) }, true);
        } else {
            writer.writeOpenTag(
                    TreeLikelihood.TREE_LIKELIHOOD,
                    new Attribute[] { new Attribute.Default("id", "treeLikelihood") }
            );
            writer.writeTag(SitePatternsParser.PATTERNS, new Attribute[] { new Attribute.Default("idref", "patterns") }, true);
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute[] { new Attribute.Default("idref", "treeModel") }, true);
            writer.writeTag(GammaSiteModel.SITE_MODEL, new Attribute[] { new Attribute.Default("idref", "siteModel") }, true);
        }
        if (clockModel == STRICT_CLOCK) {
            writer.writeTag(StrictClockBranchRates.STRICT_CLOCK_BRANCH_RATES, new Attribute[] { new Attribute.Default("idref", "branchRates") }, true);
        } else {
            writer.writeTag(DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES, new Attribute[] { new Attribute.Default("idref", "branchRates") }, true);
        }

        writer.writeCloseTag(TreeLikelihood.TREE_LIKELIHOOD);
    }

    /**
     * Generate tmrca statistics
     * @param writer the writer
     */
    public void writeTMRCAStatistics(XMLWriter writer) {

        writer.writeText("");
        for (int i = 0; i < taxonSets.size(); i++) {
            TaxonList taxa = (TaxonList)taxonSets.get(i);
            writer.writeOpenTag(
                    "tmrcaStatistic",
                    new Attribute[] {
                            new Attribute.Default("id", "tmrca(" + taxa.getId() + ")"),
                    }
            );
            writer.writeOpenTag("mrca");
            writer.writeTag("taxa", new Attribute[] { new Attribute.Default("idref", taxa.getId()) }, true);
            writer.writeCloseTag("mrca");
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute[] { new Attribute.Default("idref", "treeModel") }, true);
            writer.writeCloseTag("tmrcaStatistic");
        }
    }

    /**
     * Write the operator schedule XML block.
     * @param operators the list of operators
     * @param writer the writer
     */
    public void writeOperatorSchedule(ArrayList operators, XMLWriter writer) {
        writer.writeOpenTag(
                SimpleOperatorSchedule.OPERATOR_SCHEDULE,
                new Attribute[] { new Attribute.Default("id", "operators") }
        );

        Iterator iter = operators.iterator();
        while (iter.hasNext()) {
            Operator operator = (Operator)iter.next();
            writeOperator(operator, writer);
        }

        writer.writeCloseTag(SimpleOperatorSchedule.OPERATOR_SCHEDULE);
    }

    private void writeOperator(Operator operator, XMLWriter writer) {
        if (operator.type.equals(SCALE)) {
            writeScaleOperator(operator, writer);
        } else if (operator.type.equals(RANDOM_WALK)) {
            writeRandomWalkOperator(operator, writer);
        } else if (operator.type.equals(UP_DOWN)) {
            writeUpDownOperator(operator, writer);
        } else if (operator.type.equals(SCALE_ALL)) {
            writeScaleAllOperator(operator, writer);
        } else if (operator.type.equals(CENTERED_SCALE)) {
            writeCenteredOperator(operator, writer);
        } else if (operator.type.equals(DELTA_EXCHANGE)) {
            writeDeltaOperator(operator, writer);
        } else if (operator.type.equals(INTEGER_DELTA_EXCHANGE)) {
            writeIntegerDeltaOperator(operator, writer);
        } else if (operator.type.equals(SWAP)) {
            writeSwapOperator(operator, writer);
        } else if (operator.type.equals(UNIFORM)) {
            writeUniformOperator(operator, writer);
        } else if (operator.type.equals(SUBTREE_SLIDE)) {
            writeSubtreeSlideOperator(operator, writer);
        } else if (operator.type.equals(NARROW_EXCHANGE)) {
            writeNarrowExchangeOperator(operator, writer);
        } else if (operator.type.equals(WIDE_EXCHANGE)) {
            writeWideExchangeOperator(operator, writer);
        } else if (operator.type.equals(WILSON_BALDING)) {
            writeWilsonBaldingOperator(operator, writer);
        }
    }

    private void writeScaleOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(
                ScaleOperator.SCALE_OPERATOR,
                new Attribute[] {
                        new Attribute.Default("scaleFactor", Double.toString(operator.tuning)),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight)),
                });
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", operator.parameter1.getName()), true);
        writer.writeCloseTag(ScaleOperator.SCALE_OPERATOR);
    }

    private void writeRandomWalkOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(
                "randomWalkOperator",
                new Attribute[] {
                        new Attribute.Default("windowSize", Double.toString(operator.tuning)),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight)),
                });
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", operator.parameter1.getName()), true);
        writer.writeCloseTag("randomWalkOperator");
    }

    private void writeScaleAllOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(
                ScaleOperator.SCALE_OPERATOR,
                new Attribute[] {
                        new Attribute.Default("scaleFactor", Double.toString(operator.tuning)),
                        new Attribute.Default("scaleAll", "true"),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight)),
                });
        writer.writeOpenTag(CompoundParameter.COMPOUND_PARAMETER);
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter1.getName()) }, true);
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter2.getName()) }, true);
        writer.writeCloseTag(CompoundParameter.COMPOUND_PARAMETER);
        writer.writeCloseTag(ScaleOperator.SCALE_OPERATOR);
    }

    private void writeUpDownOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(UpDownOperator.UP_DOWN_OPERATOR,
                new Attribute[] {
                        new Attribute.Default("scaleFactor", Double.toString(operator.tuning)),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight))
                }
        );

        writer.writeOpenTag(UpDownOperator.UP);
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter1.getName()) }, true);
        writer.writeCloseTag(UpDownOperator.UP);

        writer.writeOpenTag(UpDownOperator.DOWN);
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter2.getName()) }, true);
        writer.writeCloseTag(UpDownOperator.DOWN);

        writer.writeCloseTag(UpDownOperator.UP_DOWN_OPERATOR);
    }

    private void writeCenteredOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(CenteredScaleOperator.CENTERED_SCALE,
                new Attribute[] {
                        new Attribute.Default(CenteredScaleOperator.SCALE_FACTOR, Double.toString(operator.tuning)),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight))
                }
        );
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter1.getName()) }, true);
        writer.writeCloseTag(CenteredScaleOperator.CENTERED_SCALE);
    }

    private void writeDeltaOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(DeltaExchangeOperator.DELTA_EXCHANGE,
                new Attribute[] {
                        new Attribute.Default(DeltaExchangeOperator.DELTA, Double.toString(operator.tuning)),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight))
                }
        );
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter1.getName()) }, true);
        writer.writeCloseTag(DeltaExchangeOperator.DELTA_EXCHANGE);
    }

    private void writeIntegerDeltaOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(DeltaExchangeOperator.DELTA_EXCHANGE,
                new Attribute[] {
                        new Attribute.Default(DeltaExchangeOperator.DELTA, Integer.toString((int)operator.tuning)),
                        new Attribute.Default("integer", "true"),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight)),
                        new Attribute.Default("autoOptimize", "false")
                }
        );
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter1.getName()) }, true);
        writer.writeCloseTag(DeltaExchangeOperator.DELTA_EXCHANGE);
    }

    private void writeSwapOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(SwapOperator.SWAP_OPERATOR,
                new Attribute[] {
                        new Attribute.Default("size", Integer.toString((int)operator.tuning)),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight)),
                        new Attribute.Default("autoOptimize", "false")
                }
        );
        writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", operator.parameter1.getName()) }, true);
        writer.writeCloseTag(SwapOperator.SWAP_OPERATOR);
    }

    private void writeUniformOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag("uniformOperator",
                new Attribute.Default("weight", Integer.toString((int)operator.weight)));
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref", operator.parameter1.getName()), true);
        writer.writeCloseTag("uniformOperator");
    }

    private void writeNarrowExchangeOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(ExchangeOperator.NARROW_EXCHANGE,
                new Attribute.Default("weight", Integer.toString((int)operator.weight)));
        writer.writeTag(TreeModel.TREE_MODEL, new Attribute[] { new Attribute.Default("idref", "treeModel") }, true);
        writer.writeCloseTag(ExchangeOperator.NARROW_EXCHANGE);
    }

    private void writeWideExchangeOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(ExchangeOperator.WIDE_EXCHANGE,
                new Attribute.Default("weight", Integer.toString((int)operator.weight)));
        writer.writeTag(TreeModel.TREE_MODEL, new Attribute[] { new Attribute.Default("idref", "treeModel") }, true);
        writer.writeCloseTag(ExchangeOperator.WIDE_EXCHANGE);
    }

    private void writeWilsonBaldingOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(WilsonBalding.WILSON_BALDING,
                new Attribute.Default("weight", Integer.toString((int)operator.weight)));
        writer.writeTag(TreeModel.TREE_MODEL, new Attribute[] { new Attribute.Default("idref", "treeModel") }, true);
        if (nodeHeightPrior == CONSTANT) {
            writeNodeHeightPriorModelRef(writer);
        }
        writer.writeCloseTag(WilsonBalding.WILSON_BALDING);
    }

    private void writeSubtreeSlideOperator(Operator operator, XMLWriter writer) {
        writer.writeOpenTag(SubtreeSlideOperator.SUBTREE_SLIDE,
                new Attribute[] {
                        new Attribute.Default("size", Double.toString(operator.tuning)),
                        new Attribute.Default("gaussian", "true"),
                        new Attribute.Default("weight", Integer.toString((int)operator.weight))
                }
        );
        writer.writeTag(TreeModel.TREE_MODEL, new Attribute[] { new Attribute.Default("idref", "treeModel") }, true);
        writer.writeCloseTag(SubtreeSlideOperator.SUBTREE_SLIDE);
    }

    /**
     * Write the timer report block.
     * @param writer the writer
     */
    public void writeTimerReport(XMLWriter writer) {
        writer.writeOpenTag("report");
        writer.writeOpenTag("property", new Attribute.Default("name", "timer"));
        writer.writeTag("object", new Attribute.Default("idref", "mcmc"), true);
        writer.writeCloseTag("property");
        writer.writeCloseTag("report");
    }

    /**
     * Write the trace analysis block.
     * @param writer the writer
     */
    public void writeTraceAnalysis(XMLWriter writer) {
        writer.writeTag(
                "traceAnalysis",
                new Attribute[] {
                        new Attribute.Default("fileName", logFileName)
                },
                true
        );
    }

    /**
     * Write the MCMC block.
     * @param writer the writer
     */
    public void writeMCMC(XMLWriter writer) {
        writer.writeOpenTag(
                "mcmc",
                new Attribute[] {
                        new Attribute.Default("id", "mcmc"),
                        new Attribute.Default("chainLength", chainLength+""),
                        new Attribute.Default("autoOptimize", autoOptimize?"true":"false")
                });

        writer.writeOpenTag(CompoundLikelihood.POSTERIOR, new Attribute.Default("id", "posterior"));

        // write prior block
        writer.writeOpenTag(CompoundLikelihood.PRIOR, new Attribute.Default("id", "prior"));

        writeParameterPriors(writer);

        if (nodeHeightPrior == YULE || nodeHeightPrior == BIRTH_DEATH) {
            writer.writeTag(SpeciationLikelihood.SPECIATION_LIKELIHOOD, new Attribute.Default("idref", "speciation"), true);
        } else if (nodeHeightPrior == SKYLINE) {
            writer.writeTag(BayesianSkylineLikelihood.SKYLINE_LIKELIHOOD, new Attribute.Default("idref", "skyline"), true);
        } else {
            writer.writeTag(CoalescentLikelihood.COALESCENT_LIKELIHOOD, new Attribute.Default("idref","coalescent"), true);
        }

        if (nodeHeightPrior == LOGISTIC) {
            writer.writeTag(BooleanLikelihood.BOOLEAN_LIKELIHOOD, new Attribute.Default("idref", "booleanLikelihood1"), true);
        }

        if (nodeHeightPrior == SKYLINE) {
            writer.writeTag(ExponentialMarkovModel.EXPONENTIAL_MARKOV_MODEL, new Attribute.Default("idref", "eml1"), true);
        }

        writer.writeCloseTag(CompoundLikelihood.PRIOR);

        // write likelihood block
        writer.writeOpenTag(CompoundLikelihood.LIKELIHOOD, new Attribute.Default("id", "likelihood"));

        boolean nucs = alignment.getDataType() == Nucleotides.INSTANCE;
        if (nucs && codonHeteroPattern != null) {
            for (int i = 1; i <= partitionCount; i++) {
                writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref", "treeLikelihood" + i), true);
            }
        } else {
            writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref", "treeLikelihood"), true);
        }

        writer.writeCloseTag(CompoundLikelihood.LIKELIHOOD);


        writer.writeCloseTag(CompoundLikelihood.POSTERIOR);

        writer.writeTag(SimpleOperatorSchedule.OPERATOR_SCHEDULE, new Attribute.Default("idref", "operators"), true);

        // write log to screen
        writer.writeOpenTag(MCLogger.LOG,
                new Attribute[] {
                        new Attribute.Default("id", "screenLog"),
                        new Attribute.Default(MCLogger.LOG_EVERY, echoEvery+"")
                });
        writeScreenLog(writer);
        writer.writeCloseTag(MCLogger.LOG);

        // write log to file
        if (logFileName == null) {
            logFileName = fileNameStem + ".log";
        }
        writer.writeOpenTag(MCLogger.LOG,
                new Attribute[] {
                        new Attribute.Default("id", "fileLog"),
                        new Attribute.Default(MCLogger.LOG_EVERY, logEvery+""),
                        new Attribute.Default(MCLogger.FILE_NAME, logFileName)
                });
        writeLog(writer);
        writer.writeCloseTag(MCLogger.LOG);

        // write tree log to file
        if (treeFileName == null) {
            if (substTreeLog) {
                treeFileName = fileNameStem + "(time).trees";
            } else {
                treeFileName = fileNameStem + ".trees";
            }
        }
        writer.writeOpenTag(TreeLogger.LOG_TREE,
                new Attribute[] {
                        new Attribute.Default("id", "treeFileLog"),
                        new Attribute.Default(TreeLogger.LOG_EVERY, logEvery+""),
                        new Attribute.Default(TreeLogger.NEXUS_FORMAT, "true"),
                        new Attribute.Default(TreeLogger.FILE_NAME, treeFileName)
                });
        writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
        if (clockModel != STRICT_CLOCK) {
            writer.writeTag(DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES, new Attribute[] { new Attribute.Default("idref", "branchRates") }, true);
        }
        writer.writeTag("posterior", new Attribute.Default("idref", "posterior"), true);
        writer.writeCloseTag(TreeLogger.LOG_TREE);

//        if (mapTreeLog) {
//            // write tree log to file
//            if (mapTreeFileName == null) {
//                mapTreeFileName = fileNameStem + ".MAP.tree";
//            }
//            writer.writeOpenTag("logML",
//                    new Attribute[] {
//                        new Attribute.Default(TreeLogger.FILE_NAME, mapTreeFileName)
//                    });
//            writer.writeOpenTag("ml");
//            writer.writeTag(CompoundLikelihood.POSTERIOR, new Attribute.Default("idref", "posterior"), true);
//            writer.writeCloseTag("ml");
//            writer.writeOpenTag("column", new Attribute[] {
//                        new Attribute.Default("label", "MAP tree")
//                    });
//            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
//            writer.writeCloseTag("column");
//            writer.writeCloseTag("logML");
//        }

        if (substTreeLog) {
            // write tree log to file
            if (substTreeFileName == null) {
                substTreeFileName = fileNameStem + "(subst).trees";
            }
            writer.writeOpenTag(TreeLogger.LOG_TREE,
                    new Attribute[] {
                            new Attribute.Default("id", "substTreeFileLog"),
                            new Attribute.Default(TreeLogger.LOG_EVERY, logEvery+""),
                            new Attribute.Default(TreeLogger.NEXUS_FORMAT, "true"),
                            new Attribute.Default(TreeLogger.FILE_NAME, substTreeFileName),
                            new Attribute.Default(TreeLogger.BRANCH_LENGTHS, TreeLogger.SUBSTITUTIONS)
                    });
            writer.writeTag(TreeModel.TREE_MODEL, new Attribute.Default("idref", "treeModel"), true);
            if (clockModel == STRICT_CLOCK) {
                writer.writeTag(StrictClockBranchRates.STRICT_CLOCK_BRANCH_RATES, new Attribute[] { new Attribute.Default("idref", "branchRates") }, true);
            } else {
                writer.writeTag(DiscretizedBranchRates.DISCRETIZED_BRANCH_RATES, new Attribute[] { new Attribute.Default("idref", "branchRates") }, true);
            }
            writer.writeCloseTag(TreeLogger.LOG_TREE);
        }

        writer.writeCloseTag("mcmc");
    }

    /**
     * Write the priors for each parameter
     * @param writer the writer
     */
    private void writeParameterPriors(XMLWriter writer) {
        ArrayList parameters = selectParameters();
        Iterator iter = parameters.iterator();
        while (iter.hasNext()) {
            Parameter parameter = (Parameter)iter.next();
            if (parameter.priorType != NONE) {
                if (parameter.priorType != UNIFORM_PRIOR || parameter.isNodeHeight) {
                    writeParameterPrior(parameter, writer);
                }
            }
        }

    }

    /**
     * Write the priors for each parameter
     * @param parameter the parameter
     * @param writer the writer
     */
    private void writeParameterPrior(Parameter parameter, XMLWriter writer) {
        switch (parameter.priorType) {
            case UNIFORM_PRIOR:
                writer.writeOpenTag(DistributionLikelihood.UNIFORM_PRIOR,
                        new Attribute[] {
                                new Attribute.Default(DistributionLikelihood.LOWER, "" + parameter.uniformLower),
                                new Attribute.Default(DistributionLikelihood.UPPER, "" + parameter.uniformUpper)
                        });
                writeParameterIdref(writer, parameter);
                writer.writeCloseTag(DistributionLikelihood.UNIFORM_PRIOR);
                break;
            case EXPONENTIAL_PRIOR:
                writer.writeOpenTag(DistributionLikelihood.EXPONENTIAL_PRIOR,
                        new Attribute[] {
                                new Attribute.Default(DistributionLikelihood.MEAN, "" + parameter.exponentialMean),
                                new Attribute.Default(DistributionLikelihood.OFFSET, "" + parameter.exponentialOffset)
                        });
                writeParameterIdref(writer, parameter);
                writer.writeCloseTag(DistributionLikelihood.EXPONENTIAL_PRIOR);
                break;
            case NORMAL_PRIOR:
                writer.writeOpenTag(DistributionLikelihood.NORMAL_PRIOR,
                        new Attribute[] {
                                new Attribute.Default(DistributionLikelihood.MEAN, "" + parameter.normalMean),
                                new Attribute.Default(DistributionLikelihood.STDEV, "" + parameter.normalStdev)
                        });
                writeParameterIdref(writer, parameter);
                writer.writeCloseTag(DistributionLikelihood.NORMAL_PRIOR);
                break;
            case LOG_NORMAL_PRIOR:
                writer.writeOpenTag(DistributionLikelihood.LOG_NORMAL_PRIOR,
                        new Attribute[] {
                                new Attribute.Default(DistributionLikelihood.MEAN, "" + parameter.logNormalMean),
                                new Attribute.Default(DistributionLikelihood.STDEV, "" + parameter.logNormalStdev),
                                new Attribute.Default(DistributionLikelihood.OFFSET, "" + parameter.logNormalOffset)
                        });
                writeParameterIdref(writer, parameter);
                writer.writeCloseTag(DistributionLikelihood.LOG_NORMAL_PRIOR);
                break;
            case GAMMA_PRIOR:
                writer.writeOpenTag(DistributionLikelihood.GAMMA_PRIOR,
                        new Attribute[] {
                                new Attribute.Default(DistributionLikelihood.SHAPE, "" + parameter.gammaAlpha),
                                new Attribute.Default(DistributionLikelihood.SCALE, "" + parameter.gammaBeta),
                                new Attribute.Default(DistributionLikelihood.OFFSET, "" + parameter.logNormalOffset)
                        });
                writeParameterIdref(writer, parameter);
                writer.writeCloseTag(DistributionLikelihood.GAMMA_PRIOR);
                break;
            case JEFFREYS_PRIOR:
                writer.writeOpenTag(JeffreysPriorLikelihood.JEFFREYS_PRIOR);
                writeParameterIdref(writer, parameter);
                writer.writeCloseTag(JeffreysPriorLikelihood.JEFFREYS_PRIOR);
                break;
            default:
                throw new IllegalArgumentException("Unknown priorType");
        }
    }

    private void writeParameterIdref(XMLWriter writer, Parameter parameter) {
        if (parameter.isStatistic) {
            writer.writeTag("statistic", new Attribute[] { new Attribute.Default("idref", parameter.getName()) }, true);
        } else {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute[] { new Attribute.Default("idref", parameter.getName()) }, true);
        }
    }

    /**
     * Write the log
     * @param writer the writer
     */
    private void writeScreenLog(XMLWriter writer) {
        if (alignment != null) {
            writer.writeOpenTag(Columns.COLUMN,
                    new Attribute[] {
                            new Attribute.Default(Columns.LABEL, "Posterior"),
                            new Attribute.Default(Columns.DECIMAL_PLACES, "4"),
                            new Attribute.Default(Columns.WIDTH, "12")
                    }
            );
            writer.writeTag(CompoundLikelihood.POSTERIOR, new Attribute.Default("idref","posterior"), true);
            writer.writeCloseTag(Columns.COLUMN);
        }

        writer.writeOpenTag(Columns.COLUMN,
                new Attribute[] {
                        new Attribute.Default(Columns.LABEL, "Root Height"),
                        new Attribute.Default(Columns.SIGNIFICANT_FIGURES, "6"),
                        new Attribute.Default(Columns.WIDTH, "12")
                }
        );
        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","treeModel.rootHeight"), true);
        writer.writeCloseTag(Columns.COLUMN);


        if (alignment != null) {
            boolean nucs = alignment.getDataType() == Nucleotides.INSTANCE;
            if (nucs && codonHeteroPattern != null) {
                if (codonHeteroPattern.equals("112")) {
                    writer.writeOpenTag(Columns.COLUMN,
                            new Attribute[] {
                                    new Attribute.Default(Columns.LABEL, "L(codon pos 1+2)"),
                                    new Attribute.Default(Columns.DECIMAL_PLACES, "4"),
                                    new Attribute.Default(Columns.WIDTH, "12")
                            }
                    );
                    writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref","treeLikelihood1"), true);
                    writer.writeCloseTag(Columns.COLUMN);
                    writer.writeOpenTag(Columns.COLUMN,
                            new Attribute[] {
                                    new Attribute.Default(Columns.LABEL, "L(codon pos 3)"),
                                    new Attribute.Default(Columns.DECIMAL_PLACES, "4"),
                                    new Attribute.Default(Columns.WIDTH, "12")
                            }
                    );
                    writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref","treeLikelihood2"), true);
                    writer.writeCloseTag(Columns.COLUMN);
                } else if (codonHeteroPattern.equals("123")) {
                    for (int i =1; i <= 3; i++) {
                        writer.writeOpenTag(Columns.COLUMN,
                                new Attribute[] {
                                        new Attribute.Default(Columns.LABEL, "L(codon pos "+i+")"),
                                        new Attribute.Default(Columns.DECIMAL_PLACES, "4"),
                                        new Attribute.Default(Columns.WIDTH, "12")
                                }
                        );
                        writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref","treeLikelihood" + i), true);
                        writer.writeCloseTag(Columns.COLUMN);
                    }
                }
            } else {
                writer.writeOpenTag(Columns.COLUMN,
                        new Attribute[] {
                                new Attribute.Default(Columns.LABEL, "L(tree)"),
                                new Attribute.Default(Columns.DECIMAL_PLACES, "4"),
                                new Attribute.Default(Columns.WIDTH, "12")
                        }
                );
                writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref","treeLikelihood"), true);
                writer.writeCloseTag(Columns.COLUMN);
            }
        }
        if (nodeHeightPrior == YULE || nodeHeightPrior == BIRTH_DEATH) {
            writer.writeOpenTag(Columns.COLUMN,
                    new Attribute[] {
                            new Attribute.Default(Columns.LABEL, "L(speciation)"),
                            new Attribute.Default(Columns.DECIMAL_PLACES, "4"),
                            new Attribute.Default(Columns.WIDTH, "12")
                    }
            );
            writer.writeTag(SpeciationLikelihood.SPECIATION_LIKELIHOOD, new Attribute.Default("idref", "speciation"), true);
        } else {
            writer.writeOpenTag(Columns.COLUMN,
                    new Attribute[] {
                            new Attribute.Default(Columns.LABEL, "L(coalecent)"),
                            new Attribute.Default(Columns.DECIMAL_PLACES, "4"),
                            new Attribute.Default(Columns.WIDTH, "12")
                    }
            );
            if (nodeHeightPrior == SKYLINE) {
                writer.writeTag(BayesianSkylineLikelihood.SKYLINE_LIKELIHOOD, new Attribute.Default("idref", "skyline"), true);
            } else {
                writer.writeTag(CoalescentLikelihood.COALESCENT_LIKELIHOOD, new Attribute.Default("idref","coalescent"), true);
            }
        }
        writer.writeCloseTag(Columns.COLUMN);

    }

    /**
     * Write the log
     * @param writer the writer
     */
    private void writeLog(XMLWriter writer) {
        if (alignment != null) {
            writer.writeTag(CompoundLikelihood.POSTERIOR, new Attribute.Default("idref","posterior"), true);
        } else {
            writer.writeTag(CompoundLikelihood.PRIOR, new Attribute.Default("idref","prior"), true);
        }
        if (alignment != null) {
            boolean nucs = alignment.getDataType() == Nucleotides.INSTANCE;
            if (nucs && partitionCount > 1) {
                for (int i =1; i <= partitionCount; i++) {
                    writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref","treeLikelihood" + i), true);
                }
            } else writer.writeTag(TreeLikelihood.TREE_LIKELIHOOD, new Attribute.Default("idref","treeLikelihood"), true);
        }
        if (nodeHeightPrior == YULE || nodeHeightPrior == BIRTH_DEATH) {
            writer.writeTag(SpeciationLikelihood.SPECIATION_LIKELIHOOD, new Attribute.Default("idref", "speciation"), true);
        } else if (nodeHeightPrior == SKYLINE) {
            writer.writeTag(BayesianSkylineLikelihood.SKYLINE_LIKELIHOOD, new Attribute.Default("idref", "skyline"), true);
        } else {
            writer.writeTag(CoalescentLikelihood.COALESCENT_LIKELIHOOD, new Attribute.Default("idref","coalescent"), true);
        }


        if (!fixedSubstitutionRate) {
            if (clockModel == STRICT_CLOCK) {
                writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","clock.rate"), true);
            } else {
                writer.writeTag(RateStatistic.RATE_STATISTIC, new Attribute.Default("idref","meanRate"), true);
            }
        }

        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","treeModel.rootHeight"), true);

        for (int i = 0; i < taxonSets.size(); i++) {
            TaxonList taxa = (TaxonList)taxonSets.get(i);
            writer.writeTag("tmrcaStatistic", new Attribute[] { new Attribute.Default("idref", "tmrca(" + taxa.getId() + ")") }, true);
        }

        if (nodeHeightPrior == CONSTANT) {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","constant.popSize"), true);
        } else if (nodeHeightPrior == EXPONENTIAL) {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","exponential.popSize"), true);
            if (parameterization == GROWTH_RATE) {
                writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","exponential.growthRate"), true);
            } else {
                writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","exponential.doublingTime"), true);
            }
        } else if (nodeHeightPrior == LOGISTIC) {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","logistic.popSize"), true);
            if (parameterization == GROWTH_RATE) {
                writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","logistic.growthRate"), true);
            } else {
                writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","logistic.doublingTime"), true);
            }
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","logistic.t50"), true);
        } else if (nodeHeightPrior == EXPANSION) {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","expansion.popSize"), true);
            if (parameterization == GROWTH_RATE) {
                writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","expansion.growthRate"), true);
            } else {
                writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","expansion.doublingTime"), true);
            }
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","expansion.ancestralProportion"), true);
        } else if (nodeHeightPrior == SKYLINE) {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","skyline.popSize"), true);
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","skyline.groupSize"), true);
        } else if (nodeHeightPrior == YULE) {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","yule.birthRate"), true);
        } else if (nodeHeightPrior == BIRTH_DEATH) {
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","birthDeath.birthRate"), true);
            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","birthDeath.deathRate"), true);
        }

        if (alignment != null) {
            boolean nucs = alignment.getDataType() == Nucleotides.INSTANCE;
            if (nucs) {
                if (partitionCount > 1) {
                    for (int i = 1; i <= partitionCount; i++) {
                        writer.writeTag(ParameterParser.PARAMETER,
                                new Attribute.Default("idref","siteModel" + i + ".mu"), true);
                    }
                }
                if (nucSubstitutionModel == HKY) {
                    if (partitionCount > 1 && unlinkedSubstitutionModel) {
                        for (int i = 1; i <= partitionCount; i++) {
                            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","hky" + i + ".kappa"), true);
                        }
                    } else {
                        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","hky.kappa"), true);
                    }
                } else if (nucSubstitutionModel == GTR) {
                    if (partitionCount > 1 && unlinkedSubstitutionModel) {
                        for (int i = 1; i <= partitionCount; i++) {
                            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr" + i + ".ac"), true);
                            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr" + i + ".ag"), true);
                            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr" + i + ".at"), true);
                            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr" + i + ".cg"), true);
                            writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr" + i + ".gt"), true);
                        }
                    } else {
                        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr.ac"), true);
                        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr.ag"), true);
                        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr.at"), true);
                        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr.cg"), true);
                        writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","gtr.gt"), true);
                    }
                }
            }

            if (gammaHetero) {
                if (partitionCount > 1 && unlinkedHeterogeneityModel) {
                    for (int i = 1; i <= partitionCount; i++) {
                        writer.writeTag(ParameterParser.PARAMETER,
                                new Attribute.Default("idref","siteModel" + i + ".alpha"), true);
                    }
                } else {
                    writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","siteModel.alpha"), true);
                }
            }

            if (invarHetero) {
                if (partitionCount > 1 && unlinkedHeterogeneityModel) {
                    for (int i = 1; i <= partitionCount; i++) {
                        writer.writeTag(ParameterParser.PARAMETER,
                                new Attribute.Default("idref","siteModel" + i + ".pInv"), true);
                    }
                } else {
                    writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","siteModel.pInv"), true);
                }
            }
        }

        if (clockModel != STRICT_CLOCK) {
            if (!fixedSubstitutionRate) {
                if (clockModel == UNCORRELATED_EXPONENTIAL) {
                    writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","uced.mean"), true);
                } else if (clockModel == UNCORRELATED_LOGNORMAL) {
                    writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","ucld.mean"), true);
                    writer.writeTag(ParameterParser.PARAMETER, new Attribute.Default("idref","ucld.stdev"), true);
                }
            }
            writer.writeTag(RateStatistic.RATE_STATISTIC, new Attribute.Default("idref","coefficientOfVariation"), true);
            writer.writeTag(RateCovarianceStatistic.RATE_COVARIANCE_STATISTIC, new Attribute.Default("idref","covariance"), true);
        }


    }

    /**
     * fix a parameter
     * @param id the id
     * @param value the value
     */
    public void fixParameter(String id, double value) {
        Parameter parameter = (Parameter)parameters.get(id);
        if (parameter == null) { throw new IllegalArgumentException("parameter with name, "+id+", is unknown");}
        parameter.isFixed = true;
        parameter.initial = value;
    }

    /**
     * write a parameter
     * @param id the id
     * @param writer the writer
     */
    public void writeParameter(String id, XMLWriter writer) {
        Parameter parameter = (Parameter)parameters.get(id);
        if (parameter == null) { throw new IllegalArgumentException("parameter with name, "+id+", is unknown");}
        if (parameter.isFixed) {
            writeParameter(id, 1, parameter.initial, Double.NaN, Double.NaN, writer);
        } else {
            if (parameter.priorType == UNIFORM_PRIOR) {
                writeParameter(id, 1, parameter.initial, parameter.uniformLower, parameter.uniformUpper, writer);
            } else {
                writeParameter(id, 1, parameter.initial, parameter.lower, parameter.upper, writer);
            }
        }
    }

    /**
     * write a parameter
     * @param id the id
     * @param dimension the dimension
     * @param writer the writer
     */
    public void writeParameter(String id, int dimension, XMLWriter writer) {
        Parameter parameter = (Parameter)parameters.get(id);
        if (parameter == null) { throw new IllegalArgumentException("parameter with name, "+id+", is unknown");}
        if (parameter.isFixed) {
            writeParameter(id, dimension, parameter.initial, Double.NaN, Double.NaN, writer);
        } else if (parameter.priorType == UNIFORM_PRIOR) {
            writeParameter(id, dimension, parameter.initial, parameter.uniformLower, parameter.uniformUpper, writer);
        } else {
            writeParameter(id, dimension, parameter.initial, parameter.lower, parameter.upper, writer);
        }
    }

    /**
     * write a parameter
     * @param id the id
     * @param dimension the dimension
     * @param value the value
     * @param lower the lower bound
     * @param upper the upper bound
     * @param writer the writer
     */
    public void writeParameter(String id, int dimension, double value, double lower, double upper, XMLWriter writer) {
        ArrayList attributes = new ArrayList();
        attributes.add(new Attribute.Default("id", id));
        if (dimension > 1) {
            attributes.add(new Attribute.Default("dimension", dimension + ""));
        }
        if (!Double.isNaN(value)) {
            attributes.add(new Attribute.Default("value", value + ""));
        }
        if (!Double.isNaN(lower)) {
            attributes.add(new Attribute.Default("lower", lower + ""));
        }
        if (!Double.isNaN(upper)) {
            attributes.add(new Attribute.Default("upper", upper + ""));
        }

        Attribute[] attrArray = new Attribute[attributes.size()];
        for (int i =0; i < attrArray.length; i++) {
            attrArray[i] = (Attribute)attributes.get(i);
        }

        writer.writeTag(ParameterParser.PARAMETER, attrArray, true);
    }

    /**
     * Generate XML for the starting tree
     * @param writer the writer
     */
    public void writeStartingTree(XMLWriter writer) {
        if (userTree) {
            writeUserTree(tree, writer);
        } else if (upgmaStartingTree) {
            // generate a upgma starting tree
            writer.writeComment("Construct a rough-and-ready UPGMA tree as an starting tree");
            Parameter rootHeight = getParameter("treeModel.rootHeight");
            if (rootHeight.priorType != NONE) {
                writer.writeOpenTag(
                        UPGMATreeParser.UPGMA_TREE,
                        new Attribute[] {
                                new Attribute.Default("id", "startingTree"),
                                new Attribute.Default(UPGMATreeParser.ROOT_HEIGHT, "" + rootHeight.initial)
                        }
                );
            } else {
                writer.writeOpenTag(
                        UPGMATreeParser.UPGMA_TREE,
                        new Attribute[] {
                                new Attribute.Default("id", "startingTree")
                        }
                );
            }
            writer.writeOpenTag(
                    DistanceMatrixParser.DISTANCE_MATRIX,
                    new Attribute[] {
                            new Attribute.Default(DistanceMatrixParser.CORRECTION, "JC")
                    }
            );
            writer.writeOpenTag(SitePatternsParser.PATTERNS);
            writer.writeTag(AlignmentParser.ALIGNMENT, new Attribute[] { new Attribute.Default("idref", "alignment") }, true);
            writer.writeCloseTag(SitePatternsParser.PATTERNS);
            writer.writeCloseTag(DistanceMatrixParser.DISTANCE_MATRIX);
            writer.writeCloseTag(UPGMATreeParser.UPGMA_TREE);
        } else {
            // generate a coalescent tree
            writer.writeComment("Generate a random starting tree under the coalescent process");
            Parameter rootHeight = getParameter("treeModel.rootHeight");
            if (rootHeight.priorType != NONE) {
                writer.writeOpenTag(
                        CoalescentSimulator.COALESCENT_TREE,
                        new Attribute[] {
                                new Attribute.Default("id", "startingTree"),
                                new Attribute.Default("rootHeight", "" + rootHeight.initial)
                        }
                );
            } else {
                writer.writeOpenTag(
                        CoalescentSimulator.COALESCENT_TREE,
                        new Attribute[] {
                                new Attribute.Default("id", "startingTree")
                        }
                );
            }
            writer.writeTag(TaxaParser.TAXA, new Attribute[] { new Attribute.Default("idref", "taxa") }, true);
            writeInitialDemoModelRef(writer);
            writer.writeCloseTag(CoalescentSimulator.COALESCENT_TREE);
        }
    }

    public void writeInitialDemoModelRef(XMLWriter writer) {
        if (nodeHeightPrior == CONSTANT) {
            writer.writeTag(ConstantPopulationModel.CONSTANT_POPULATION_MODEL, new Attribute[] { new Attribute.Default("idref", "constant") }, true);
        } else if (nodeHeightPrior == EXPONENTIAL) {
            writer.writeTag(ExponentialGrowthModel.EXPONENTIAL_GROWTH_MODEL, new Attribute[] { new Attribute.Default("idref", "exponential") }, true);
        } else {
            writer.writeTag(ConstantPopulationModel.CONSTANT_POPULATION_MODEL, new Attribute[] { new Attribute.Default("idref", "initialDemo") }, true);
        }
    }

    public void writeNodeHeightPriorModelRef(XMLWriter writer) {
        if (nodeHeightPrior == CONSTANT) {
            writer.writeTag(ConstantPopulationModel.CONSTANT_POPULATION_MODEL, new Attribute[] { new Attribute.Default("idref", "constant") }, true);
        } else if (nodeHeightPrior == EXPONENTIAL) {
            writer.writeTag(ExponentialGrowthModel.EXPONENTIAL_GROWTH_MODEL, new Attribute[] { new Attribute.Default("idref", "exponential") }, true);
        } else if (nodeHeightPrior == LOGISTIC) {
            writer.writeTag(LogisticGrowthModel.LOGISTIC_GROWTH_MODEL, new Attribute[] { new Attribute.Default("idref", "logistic") }, true);
        } else if (nodeHeightPrior == EXPANSION) {
            writer.writeTag(ExpansionModel.EXPANSION_MODEL, new Attribute[] { new Attribute.Default("idref", "expansion") }, true);
        } else if (nodeHeightPrior == SKYLINE) {
            writer.writeTag(BayesianSkylineLikelihood.SKYLINE_LIKELIHOOD, new Attribute[] { new Attribute.Default("idref", "skyline") }, true);
        } else if (nodeHeightPrior == YULE) {
            writer.writeTag(YuleModel.YULE_MODEL, new Attribute[] { new Attribute.Default("idref", "yule") }, true);
        } else if (nodeHeightPrior == BIRTH_DEATH) {
            writer.writeTag(BirthDeathModel.BIRTH_DEATH_MODEL, new Attribute[] { new Attribute.Default("idref", "birthDeath") }, true);
        } else {
            throw new RuntimeException("No coalescent model has been specified so cannot refer to it");
        }
    }

    /**
     * Generate XML for the user tree
     * @param tree the user tree
     * @param writer the writer
     */
    private void writeUserTree(Tree tree, XMLWriter writer) {

        writer.writeComment("The starting tree.");
        writer.writeOpenTag(
                "tree",
                new Attribute[] {
                        new Attribute.Default("height", "startingTree"),
                        new Attribute.Default("usingDates", (maximumTipHeight > 0 ? "true": "false"))
                }
        );
        writeNode(tree, tree.getRoot(), writer);
        writer.writeCloseTag("tree");
    }

    /**
     * Generate XML for the node of a user tree.
     * @param tree the user tree
     * @param node the current node
     * @param writer the writer
     */
    private void writeNode(Tree tree, NodeRef node, XMLWriter writer) {

        writer.writeOpenTag(
                "node",
                new Attribute[] { new Attribute.Default("height", ""+tree.getNodeHeight(node)) }
        );

        if (tree.getChildCount(node) == 0) {
            writer.writeTag("taxon", new Attribute[] { new Attribute.Default("idref", tree.getNodeTaxon(node).getId()) }, true);
        }
        for (int i =0; i < tree.getChildCount(node); i++) {
            writeNode(tree, tree.getChild(node, i), writer);
        }
        writer.writeCloseTag("node");
    }
}

